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We investigate the structure and physical properties of the undoped bilayer silicene (BLS). Our first principle 
(FP) calculations reveal that the band structure of the system with D 3 d symmetry is intrinsically metallic with 
pocket Fermi surfaces around each K point, which opens the door to formation of a superconducting state. Our 
further study based on random phase approximation (RPA) identifies the system as a chiral d + id' superconduc- 
tor, mediated by the strong spin fluctuation on border of the antiferromagnetic spin density wave (SDW) order. 
Tunable Fermi pocket area via strain enables high superconducting critical temperature, which emerges when 
the SDW critical interaction strength is tuned near that of the real interaction. Our discovery not only sheds 
light upon the realization of the long-sought chiral d + id' superconductivity (SC), but will also bring the exotic 
unconventional SC into the familiar Si-industry. 



Chiral SC is a special kind of topological SC characterized 
by time reversal symmetry breaking 1 . In the past few years, a 
surge of theoretical proposals have been raised on the experi- 
mental realization of this kind of unconventional SC, includ- 
ing such examples as the triplet p x ± ip y (p + ip') pairing 2-4 
and the singlet d x 2_ y 2 ± id xy (d + id') pairing 5-13 . While 
the former has probably been realized in the Sr2Ru04 sys- 
tem 14 , no experimental evidence has now been detected for 
the latter although the cuprates family 5-7 and more probably 
the doped graphene 8-13 have been proposed as possible candi- 
dates. Here, we predict the realization of the singlet d+id! SC 
in the undoped BLS, which is easier than to dope the graphene 
system. As a result of its nontrivial topological property, this 
intriguing pairing state will bring a series of interesting exper- 
imental consequences such as quantized boundary current 5 , 
spontaneous magnetization 5,7 and quantized spin and thermal 
Hall conductance 7 . 

Silicene, considered as the silicon-based counterpart of 
the graphene, has attracted much attention both theoretically 
and experimentally 15-27 . Similar honeycomb lattice structure 
of the two systems let them share most of their marvellous 
physical properties, including gapless Dirac fermions at the 
Brillouin-Zone corner, and the quantum spin Hall effect when 
spin-orbital coupling turns on 16,17 . Remarkable difference be- 
tween the graphene and the silicene systems mainly lies in the 
nonplanar low-buckled (LB) lattice structure of the latter 16,18 , 
which originates from the weakened it— bond between neigh- 
boring silicon atoms. Just like the bilayer graphene (BLG), 
silicene can also take the form of its bilayer version, which has 
recently been synthesized 24 . However, due to the LB struc- 
ture of each silicene layer, there are actually different stacking 
ways between the silicene bilayer, including the one studied in 
the Ref 28 , which might not be energetically optimized. Thus, 
it is necessary to identify the realized stacking ways and the 
stability of the bilayer structure of the system. 

In this paper, we first identify the AB-bt bilayer struc- 
ture (introduced below) as the energetically most favored one 
through FP calculations, and find that the corresponding elec- 
tronic band structure is intrinsically metallic for the undoped 
case with much larger Fermi pockets than the BLG 29 , opening 
the door to formation of a superconducting state. We further 
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FIG. 1 . Geometry and phonon spectrum of the BLS. a, Optimized 
geometry of the BLS. b, The corresponding phonon spectrum. In a, 
both the top view (upper) and side view (lower) are shown. The verti- 
cal bond length l v , the nearest-neighbor in plane bond length l n , and 
the angle between them are marked, together with the hopping in- 
tegrals between each two of the four atoms Ai, Bi (i = 1, 2) within 
a unit cell. 



did a RPA based study of the system when realistic Hubbard- 
interaction is added into our tight-binding (TB) model. Sub- 
sequently, a collinear antiferromagnetic SDW order is found 
at realizable interaction parameters. Interestingly, below but 
near the SDW critical point which is tunable via strain, a chi- 
ral d + id' pairing state emerges with possibly high super- 
conducting critical temperature, mediated by the strong spin 
fluctuation on border of the SDW order. The exotic chiral 
d + id' SC in the BLS can thus be manipulated via strain, 
which opens prospects for both studying the unconventional 
topological SC in new playground and for applications in the 
Si based electronics. 

Crystal and Band Structure: The crystal and elec- 
tronic band structures of the BLS reported below are ob- 
tained through our FP calculations based on Density func- 
tional theory (DFT). The electronic band structure of the sys- 
tem is obtained self-consistently by using the projector aug- 
mented wave (PAW) pseudopotential method implemented in 
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TABLE I. TB parameters fitted for the BLS in comparison with 
those of the BLG 33 . The unit is eV. 



FIG. 2. Band structure of the BLS. a, The band structure of BLS 
corresponding to the optimized structure, b, The FS patches around 
each K point. In a, zooming in the low energy band (right) is also 
shown. In b, the central pocket (red) is electron-like and the outer 
three identical pockets (green) are hole-like, where the total areas of 
the two kinds of pockets are equal. 



the VASP package 30 . The exchange-correlation potential is 
treated by Perdew-Burke-Ernzerhof (PBE) potential 31 . 

As a consequence of the LB structure of each silicene layer, 
there are actually four stacking ways (see Supplementary in- 
formation I ) between the upper and lower layers in the sys- 
tem. Our FP calculations revealed that two of them are stable, 
among which the energetically favored one (named as the AB- 
bt structure) is shown in Fig. la, and the corresponding phonon 
spectrum 32 is shown in Fig. lb. 

From Fig. la, the bottom (Ai-sublattice) of the upper sil- 
icene layer couples with the top (^-sublattice) of the lower 
layer vertically with a bond-length l v = 2. 53 A, while the two 
sublattices (A and B) within a layer couples with a bond- 
length l n = 2.32A. Approximately equal bond lengths, to- 
gether with the bond- angle = 106.60° between the two 
bonds describes an orbital hybridization more like the sp 3 
type (with bond angle 6 = 109.47°) than the planar sp 2 type. 
From Fig. lb, the phonon frequencies obtained are real at all 
momenta, which suggests a stable structure. The energy of 
this configuration is -19.65eV per unit cell, lower than that of 
the configuration studied in the literature 28 , which is -19.51 eV 
per unit cell. It's noting here that the symmetry group of the 
system is D 3 d- 

The band structure of the BLS with AB-bt stacking way is 
shown in Fig. 2a (left), together with its low energy zooming 
in (right). The most obvious feature of this band structure is 
the 300meV overlap between the valence band and the con- 
duction band, much larger than the 1.6meV in the BLG and 
the 40meV in the graphite. Another important feature is the 
asymmetric band crossings present not only at the K point, but 
also on the K-r axis with an energy difference between them. 
Such an asymmetric band crossing results in a three-folded 
symmetric pocket Fermi surface (FS) structure surrounding 
each K point, as shown in Fig. 2b, where the central electron 
pocket is accompanied by three identical outer hole pockets 
with equal total area. Here, only the FS patches around one K 
point are present. The other patches can be obtained by six- 
folded rotations around the T point, as required by the D^d 
symmetry and the time reversal invariant of the system. 
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model in the basis ^2), \Ai) : I ^2)}, which well cap- 

tures all the low energy features of the above band structure 
near the FS, 
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Here A^Bi (i = 1,2) represent the basis mainly composed 
of the 3p 2 -orbitals localized around each of the four silicon 
atoms within a unit cell. The hopping integrals t n , t, £2 and 
ts between each two orbitals are marked in Fig. la. The phase 
factor /(k) = £ a e ik ' R «, with H a (a = 1,2,3) to be the 
nearest-neighbor vector. Finally, notice the small effective on 
site energy difference A between the A and B atoms. The fit- 
ted parameters of the system in comparison with those of the 
BLG are listed in Table I, from which the most obvious fea- 
ture of the BLS lies in the dominating role of the vertical inter- 
layer hopping t. The resulting strong bonding-antibonding en- 
ergy split between A\ and A2 orbitals pushes them far away 
from the Fermi level, leaving B\ and B2 orbitals to form a 
low energy subspace which takes responsibility for the main 
physics of the system. 

It's important to point out here that the low energy band 
structure of the system is considerably sensitive to the biaxial 
strain exerted on each silicene layer. As shown in Fig. 3c, for 
small strains which keep the symmetry and FS topology of 
the system, the total area of the electron or hole pockets feels 
a considerable variation. This tunable property of the band 
structure turns out to be very important for our following dis- 
cussions. 

RPA and SDW: Let's consider the following 4-band 
Hubbard-model of the system, 

H= ^2 c La^«/5( k ) c k/3a + U ^ n ia ^n ia ^ (2) 

kcr,Q!/3 i,oc=l,4 

where H(k) is defined by (1), i and a (f3) denote the unit cell 
and orbital indices respectively. Standard multi-orbital ran- 
dom phase approximation (RPA) 34-39 (see also Supplemen- 
tary information II ) approach is adopted in our study. 
The free susceptibility (U — 0) of the model is, 



v (0)/l,/2 



To proceed, we construct the following effective 4-band TB where U (i 



(q,r)= (^(ki^c^ki+q,^ 

ki,k 2 

c+(k 2 + q,0)Q 4 (k 2 ,0)) , (3) 
= 1, • • • ,4) denotes orbital index. The k de- 
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FIG. 3. Free static susceptibility, the SDW ordered spin pattern 
and the biaxial strain-dependence of Fermi pocket area ratio and 
the SDW critical interaction strength, a, k-dependence of the free 
static susceptibility, b, The SDW ordered spin pattern, c, The biaxial 
strain-dependence of Fermi pocket area ratio, viz, the ratio of the 
total area of the electron and hole pockets against the total area of 
Brillouin Zone, and the SDW critical value U c of the BLS. 



pendent static susceptibility of the system defined by the 
largest eigenvalue of the susceptibility matrix x[°m (q) = 

Xm]rn (q, iv — 0) is shown in Fig. 3a, which displays a dis- 
tribution centering around the r -point. 

When interaction turns on, the spin (x^) or charge (x^) 
susceptibility in the RPA level is given by 



>(c)) 



(q, iv) =\lT X (0) (q, iv) (U)] X (0) (q, iv) • (4) 



Here (U) is 16 x 16 matrix, whose only four nonzero ele- 
ments are (U)w = U (/i = 1,...,4). It's clear that the 
repulsive Hubbard-interaction suppresses x^ and enhances 
X^. When the interaction strength U is enhanced to a crit- 
ical value U c , the spin susceptibility of the model diverges, 
which implies the instability toward long-range SDW order. 
The ordered spin structure of this bilayer system determined 
by the eigenvector of the spin susceptibility matrix xfm (q) = 

Xm]ln (q, iv = 0) corresponding to its largest eigenvalue is 
shown in Fig. 3b, from which one finds an antiferromagnetic 
state with antiparallelly aligned spin patterns within a unit 
cell. The ordered moments are mainly distributed on the Bi 
(i = 1,2) atoms which take responsibility for the low energy 
physics of the system. It's noting here that with the enhance- 
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FIG. 4. k-dependence of the gap function of the d x i_ y i sym- 
metry and the interaction strength U dependence of the largest 
eigenvalue A of the linearized gap function, a, k-dependence of 
the gap function of the d x i_ y i symmetry, which is one of the lead- 
ing symmetry of the system (the other one, i.e. the d xy is shown 
in the Supplementary information III), b, The interaction strength U 
dependence of the largest eigenvalue A of the linearized gap function 
(6), which is related to T c through T c = l.l3hw D e~ 1/x . In a, the 
gap function is antisymmetric about the axis x — ±y shown in the 
reciprocal space. In b, results for different strain values are shown. 



ment of the strain and hence the Fermi pocket area, the SDW 
critical value U c feels an obvious variation from the 1 .48eV at 
zero strain to the 1.1 8eV at the strain of 0.06. Such a range is 
probably realizable for the Hubbard- U of the 3p z orbitals of 
the silicon atoms. 

Superconductivity: Through exchanging antiferromag- 
netic spin fluctuations between each cooper pair, unconven- 
tional chiral d + id! SC emerges in the BLS system. 

Let's consider the process when a cooper pair with mo- 
menta (k 7 , — k') in the /3-th (/? = 1,2) band are scattered to 
(k, — k) in the a-th (a = 1, 2) band by charge or spin fluctu- 
ations. This process results in the following effective interac- 
tion, 



eff 



V^(k, 104 (k) C t (-k) C 0(-k%(k'), (5) 
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where the effective interaction vertex V^(k, k 7 ) determined 
by the RPA susceptibility (4) can be found in Ref 39 (see also 
Supplementary information II). 

From this effective interaction, we obtained the following 
linearized gap equation 39 near T c , 

Here, the integration is along various FS patches labelled by a 
or p. The i^(k') is Fermi velocity at k' on the f3— th FS patch, 
and k'u represents the component along that patch. In this 
eigenvalue problem, the normalized eigenvector A a (k) rep- 
resents the relative value of the gap function on the a—th FS 
patches near T c which is related to the eigenvalue A through 
T c = 1.13fkJc>e _1 / A . Here Hljd is a typical energy scale for 
the spin or charge fluctuation, approximated as the low en- 
ergy band width, i.e. hup « 300meV. The leading pairing 
symmetry of the system is thus determined by the eigenvector 
A a (k) corresponding to the largest eigenvalue. 

Our RPA calculations on the BLS identify an exactly de- 
generate d xy and d x i_ y i doublets as the leading pairing sym- 
metry of the system for U < U c at all strain values, which is 
robust against small doping (see Supplementary information 
III). Both symmetries are singlet with nodal gap functions. 
While the d x i_ y i shown in Fig.4a is antisymmetric about the 
axis x = ±y in the reciprocal space, the d xy shown in Sup- 
plementary information III is symmetric about them. The two 
gap functions form a 2D E g representation of the D 3 d point- 
group of the system. For both symmetries, two gap nodes 
are present on each Fermi pocket. The U dependence of the 
eigenvalue A which is related to T c is shown in Fig.4b for 
different strains. Clearly, T c increases with the Hubbard-U 
and rises promptly at U/U c < 1 as a result of the strongly 
enhanced antiferromagnetic spin fluctuation near the critical 
point. Since U c is tunable via strain as shown in Fig. 3c, the 
ratio U/U c varies within a range which provides basis for the 
realization of the relation U/U c < 1 which is crucial for the 
high T c of the system. For example, for A « 0.3 attainable 
by different strains shown in Fig.4b, the T c obtained can be as 
high as over 80K, although it is usually overestimated in the 
RPA level. For real material, whether high T c can be acquired 
is determined by how near U/U c can be tuned to 1. 

Since the two d-wave pairing symmetries are degener- 
ate, one would conjecture that they would probably be 
superposed 11 to lower the energy below T c . To determine how 
they are superposed at the ground state, let's set the gap func- 
tion as A (k) = Ai Ad w2 _ v 2 ( k ) + A 2 A dxy (k), where the 
superposition coefficients A* (i = 1,2) are determined by 
minimization of the mean- field ground state energy of the ef- 
fective Hamiltonian H eff = H band + V eff . Here H band is 
the kinetic part of (2) and V e ff is given by (5). Our energy 
minimization gives A 2 = zAi, which just leads to the long 
sought nodeless chiral d + id' SC! This superposition manner 
between the two d-wave pairings satisfies the requirement that 
the gap nodes should avoid the FS to lower the energy. With 
intrinsically complex gap functions, this pairing breaks time 



reversal symmetry and will bring a lot of exotic properties. It 
is just a singlet analogy of the extensively studied p + ip' SC. 

Our RPA calculations for the system also identify a possi- 
ble nodeless f-wave pairing to be the leading symmetry in the 
triplet channel. This pairing also breaks time reversal symme- 
try and the gap functions change sign with every 60° rotation, 
which belongs to A\ u representation of D 3 d (see Supplemen- 
tary information III) . However, its T c is much lower than that 
of the d + id f pairing. 

Conclusion: We have performed a FP calculation on the 
BLS. Through energy optimization, we identified an D 3 d 
symmetric AB-bt stacking structure for the BLS. The band 
structure corresponding to this crystal structure is intrinsically 
metallic, with Fermi pockets around each K point whose areas 
are tunable via strain. We have further carried out a system- 
atic RPA based study for the system. Our study suggests that 
above realizable interaction strength, collinear antiferromag- 
netic order will develop in the system. Below the SDW critical 
point U c , a chiral d + id' SC emerges in the system, whose T c 
rises to its maximum at U < U c as a result of the strongly 
enhanced antiferromagnetic spin fluctuations near the critical 
point. Tunable U c with strain makes high temperature SC pos- 
sible in this system, the realization of which will bring a new 
epoch to the familiar Si-industry. 
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TABLE A.l. Optimized energies and geometric parameters cor- 
responding to the four possible structures of the BLS. The second 
column gives the energy per unit cell. In the third column, the first 
and second number provide the bond-lengths of the nearest neighbor 
intra-layer and inter-layer valence bonds respectively. The final col- 
umn is the angle between these two kinds of valence bonds. Note 
that in the AB-bb structure, the exchange symmetry between the up- 
per and lower layers has been broken, as reflected in the different 
bond-lengths and bond-angles for the two layers. 



Structure Energy(eV) Bond-length(A) Bond-Angle 



AA-bb 
AA-bt 
AB-bt 
AB-bb 



-19.14 
-19.51 
-19.65 
-19.30 



2.28,5.12 
2.32, 2.46 
2.32, 2.53 
2.31(2.34), 2.92 



101.52° 
106.48° 
106.60° 

107.77°(109.53°) 



A. POSSIBLE STRUCTURES 
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Due to the weakened tt— bond between adjacent Si-atoms, 
the silicene layer possesses a low-buckled non-planar struc- 
ture, which is different from that of the graphene. In this struc- 
ture, a silicene layer is divided into two halves, i.e. the top 
and bottom, with each occupying one of the two sublattices 
(A and B) of the original honeycomb lattice. As a result of 
this non-planar structure of a silicene layer, there can be four 
possible stacking ways between the two layers in the bilayer 
silicene (BLS), i.e. the AA bottom-bottom (AA-bb), the AA 
bottom-top (AA-bt), the AB bottom-top (AB-bt), and the AB 
bottom-bottom (AB-bb). The geometries and phonon spectra 
corresponding to these four structures are shown in Fig.A.l, 
and the corresponding optimized energies and parameters are 
listed in Table A.l 

From Fig.A.l, it's clear that the four possible structures 
of the system possess different symmetries. The AA-bb and 
AB-bt structures shown in Fig.A.la and Fig.A.lc possess the 
Dsd symmetry, while the point group of the AA-bt structure 
shown in Fig. A. lb is D 3 h- The AB-bb structure shown in 
Fig. A. Id has the smallest symmetry group among the four, 
i.e. Cs v , since the different bond-lengths and bond-angles be- 
tween the two layers shown in Fig.A.le and Table A.l. From 
the phonon spectra corresponding to these structures, one 
finds that Fig.A.lb and Fig.A.ld are unstable toward struc- 
ture reconfiguration because of the imaginary frequencies ob- 
tained, leaving only the two symmetric structures shown in 
Fig. A. la and Fig. A. lc as possible candidates of realized struc- 
ture. 

The AA-bb structure shown in Fig.A.la can be considered 
as a simple analogy to the AA-stacked BLG which consists 
of two graphene layers weakly coupled through the van der 
Waals interaction. Consequently, one finds from Table A.l 
that the inter-layer bond-length corresponding to this struc- 
ture is extremely large. This structure is energetically worst 
although it is locally stable. The AA-bt structure shown in 
Fig. A. lb is just the one studied in Ref 28 . Although the energy 
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FIG. A. 1 . Geometries and phonon spectra corresponding to the 
four possible structures of the BLS. a-d, Geometries and phonon 
spectra for the AA-bb (a), AA-bt (b), AB-bt (c) and AB-bb (d) struc- 
tures of the BLS. In each show, the left up is the top view, the left 
down is the side view and the right is the phonon spectrum corre- 
sponding to each structure. Note that negative frequencies shown 
actually represent imaginary values obtained and thus imply struc- 
ture instability, e, More details for the configuration of the AB-bb 
structure, where one finds the breaking of the exchange symmetry 
between the two layers in this structure. 
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TABLE A.2. TB parameters fitted for the BLS via different 
strain. The unit is eV. 



Strain 
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t 2 


ts 
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0.00 


1.130 
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0.158 
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-0.100 



of this structure is somewhat favorable, it is locally unstable. 
The structure Fig. A. Id is not energetically favored either. The 
realized structure of the BLS should be the AB-bt structure 
shown in Fig.A.lc which is focused in the present work as it 
is not only energetically most favored but also locally stable. 
The bond- angle of this structure is 106.60°, very near to that 
of the sp 3 hybridization, i.e 109.47°. The inter-layer bond- 
length is 2. 53 A, near to the 2. 32 A of the intra-layer bond- 
length. Combination of the two features describes an sp 3 -like 
orbital hybridization rather than a planar sp 2 -like one. 

The band structures corresponding to the above four crystal 
structures are shown in Fig. A. 2. In the AA-bb band struc- 
ture shown in Fig.A.2a, two crossings slightly deviating from 
each K point are present at the Fermi level, which leads to 
a semi-metal ground state. This state is not favored by en- 
ergy, as illustrated above. The AA-bt band structure shown in 
Fig. A. 2b has been studied in Ref 28 . The two crossings locat- 
ing on the F— M and T— K lines are present at the Fermi level, 

I 



J 



where a,j3 = 1, 4 are band indices, eg and (k) are the 
a— th eigenvalue and eigenvector of the H{k) matrix (equa- 
tion (1) in the main text) respectively and rip is the Fermi- 



also leading to a semi-metal state. However, our FP calcula- 
tions revealed that this state is neither energetically most fa- 
vored nor stable. The AB-bt band structure shown in Fig.A.2c 
is the finally realized one. Large band overlap and asymmetric 
band crossings are present in this band structure, which leads 
to the pocket Fermi surface (FS) structure shown in Fig.A.2e. 
The six-folded rotation symmetry of the FS is required by the 
Dsd and time reversal symmetry of the system. From this 
realized band structure, one obtains an intrinsically metallic 
ground state with finite DOS at the Fermi level, which opens 
the door to formation of SC. The AB-bb band structure shown 
in Fig.A.2d is somewhat similar with Fig.A.2c in that it also 
possesses band overlap. However, it is energetically unfa- 
vored. 

It's important that the biaxial strain exerted on each silicene 
layer has obvious influence on the low energy band struc- 
ture of the system. In the following Table A. 2, we list the 
tight-binding (TB) parameters fitted for each strain up to 0.06, 
under which the symmetry and FS topology of the system is 
kept. 

B. FOUR-ORBITAL RPA 

In this part, we provide details of the four-orbital ran- 
dom phase approximation (RPA) approach 34-39 adopted in the 
study of the above Hubbard-model (formula (2) in the main 
text). 



A. RPA susceptibility 

Let's define, 



(B.l) 



I 

to be the free susceptibility for U = 0, with ^ (i = 1, • • • , 4) 
denoting orbital index. The explicit formulism of x^ in the 
momentum-frequency space is, 



(B.2) 



I 

Dirac distribution function. 

When interaction turns on, we define the spin (X (s) ) and 
charge (x^) susceptibility as follow, 



xl 0) ^ 2 (q,r)^ £ (^(^,^(^+^^+(^ + ^0)0^(^,0))^ 

ki,k 2 



y (0)ii 



(q, iu n ) 



± e gwcortg a* + <off-(k + q) n : (£k : q i" n f v 



k,a,/3 



S k+q 



7 



(c)d 


h (q,r) 


(s z )h 

x h,h 


h (q,r) 


(s+-)h 


/2 (q,r) 


(s-+)h 


/2 (q,r) 



£'' 2 (q^) = ^ E (^^^(k^r)^,.^^ +q,r)C+ ff2 (k 2 + q,0)C ;4 , ff2 (k 2 ,0) 

ki ,k2 ,<Ji ,<J2 

"i ,1,,a (q,T) = J E ^ (^^(ki.rJC.^Ckx + q, r)C+ (Ta (k a + q,0)C l4 , (T2 (k a) 0)), 



ki,k 2 ,cr 1 ,cr 2 



ki,k 2 



ki,k 2 



(B.3) 



Note that in non-magnetic state we have = ^ = 

X^ s + ) = X^> an d when J7 = we have x^ = X^ = X^ • 
In the RPA level, the spin/charge susceptibility for the 
Hubbard-model (formula (2) in the main text) is, 

X (s) (q, iv) =[l- X (0) (q, %v) (U)] X (0) (q, «/) , 

X (c) (q, zi/) = [/ + X (0) (q, ^) (tf)] _1 X (0) (q, ^) (B.4) 

where x^ c ^ (q, zV n ), X^(q?^ n ) and (Z7) are 16 x 16 
matrices with elements of the matrix (U) to be (U)q^ = 
115^=^=0=^. 



B. Pairing symmetry 



Let's consider a cooper pair with momentum/orbital 
(k% — k's), which could be scattered to (kp, — kg) by charge 
or spin fluctuations. In the RPA level, The effective interac- 
tion induced by this process is as follow, 

V e R ff A = E C(fc,fc')ct(k) C t(-k) Cs (-k')c t (k'), 

(B.5) 

In the singlet channel, the effective vertex I^f (fc, A: 7 ) is given 
as follow, 



1 st 



(fc, k') = (U)* + \{(U) [3 X (S) (k - k') - X (c) (k - k')] (U)}* + 
\ {(U) [3 X (S) (k + k') - x {c) (k + k')] (U)}^ , 



while in the triplet channel, it is 



1 st 



(k, k') = -\ {([/) [ X (s) (k - k') + x (c) (fc - k 1 )] {u)Y + 
\ {(U) [x is) (k + k') + x (c) (k + k')] (U)}^ 



(B.6) 



(B.7) 



Notice that the vertex T p s 1 (k, k') has been symmetrized for the Projecting the above effective interaction (B.5) into the two 

singlet case and anti-symmetrized for the triplet case. Gener- bands which cross the FS, we obtain the following low energy 

ally we neglect the frequency-dependence of T and replace it effective Hamiltonian for the cooper pairs near the FS, 
byr^(fc,fc')«C(k,k',0). 



v *ff= E ^(k,k')4(k)4(-k) c ^(-k')^(k') ) 

a/3,kk' 



where a, P = 1, 2 and V a(3 (k, k') is 



V^(k,k') = Re ^(k,^^)^'^^^'^-^^^-^)^^^)- 

pqst, kk' 



(B.8) 



(B.9) 
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Here only intra-band pairing is considered for small U. 

From the low energy effective Hamiltonian (B.8), one can 
obtain the following linearized gap equation 39 to determine 
the T c and the leading pairing symmetry of the system, 

(B.10) 

Here, the integration and summation is along various FS 
patches labelled by a or f3. The i^(k') is Fermi velocity at k' 
on the f3— th FS patch, and fey represents the component along 
that patch. In this eigenvalue problem, the normalized eigen- 
vector A a (k) represents the relative value of the gap function 
on the a— th FS patches near T c determined by 

T c = 1.13Huj D e- 1/x . (B.ll) 

The leading pairing symmetry of the system is thus deter- 
mined by the largest eigenvalue A of Eq.(B.lO). 



C. MORE DETAILS OF THE PAIRING SYMMETRIES 
OBTAINED 

Our RPA calculations reveal that at half-filling, the leading 
pairing symmetries in the BLS are the d x i_ y i and d xy dou- 



blets. The two gap functions are symmetry related: one is 
obtained from the other through the symmetry operations of 
Dsd, and thus they are exactly degenerate. While the normal- 
ized relative gap function of the d x i_ y i symmetry has been 
shown in Fig. 4a of the main text, here we provide that of the 
d xy symmetry in Fig. C. la. Obviously, the gap function for 
this pairing is symmetric about the axis x = ±y shown in the 
reciprocal space. 

It is interesting that our RPA calculations also identify a 
possible f-wave SC shown in Fig. C. lb to be the leading one 
in the triplet channel at half-filling. The gap function of this 
pairing changes sign with every 60° rotation. It is antisymmet- 
ric about the three diagonal lines shown in the Brillouin-Zone, 
which are just the gap node lines. This gap function belongs 
to Aiu irreducible representation of Dsd- Similarly with the 
d + id' symmetry, this pairing also satisfies the requirement 
that gap nodes should avoid the FS. Positive eigenvalue A cor- 
responding to this symmetry suggests finite T c of the pairing, 
which is nevertheless much lower than that of the d + id' sym- 
metry. 

The d + id' pairing found here is robust against small dop- 
ing. Our RPA calculations yield that for U = leV, the leading 
pairing symmetry of the system is still d + id! either for the 
5% electron-doped case shown in Fig. C. 2a or for the 5% hole- 
doped case shown in Fig. C. 2b although the FS topology has 
been modified. 
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FIG. C.l. Normalized relative gap function of the leading pairing 
symmetries in different channels for the undoped case, a, Gap 

function of the d xy symmetry, which is the other leading symmetry 
in addition to the d x 2_ y 2 pairing shown in Fig.4a of the main text, b, 
Gap function of the f-wave pairing obtained as the leading symmetry 
in the triplet channel. 




FIG. C.2. Leading pairing symmetries for the doped cases. a,b, 

Leading pairing symmetries of the 5% electron-doped case (a) and 
5% hole-doped case (b) for U = leV. In each show, the left is the 
d x 2_ y 2 and the right is the d xy . 



